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We investigate the possibility of quantum (or wave) chaos for the Bogohubov excitations of a 
Bose-Einstein condensate in bilhards. Because of the mean field interaction in the condensate, the 
Bogoliubov excitations are very different from the single particle excitations in a non-interacting 
system. Nevertheless, we predict that the statistical distribution of level spacings is unchanged by 
mapping the non-Hermitian Bogoliubov operator to a real symmetric matrix. We numerically test 
our prediction by using a phase shift method for calculating the excitation energies. 
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In recent years, the realization of Bose-Einstein con- 
densation (BEG) of dilute gases jl| has opened new op- 
portunities for studying dynamical systems in the pres- 
ence of many-body interactions. However, most previous 
investigations have focused on one dimensional or high 
dimensional separable systems and the dynamics of BEG 
in nonseparable systems with two or more degrees of free- 
dom have not received much attention . 

In the hnear Schrodinger equation, systems with two 
or more degrees of freedom can be characterized by the 
statistics of energy levels: the typical distribution of the 
spacing of neighboring levels is Poisson or Gaussian En- 
sembles for separable or nonseparable systems respec- 
tively 0. In the limit of short wavelengths (geomet- 
ric optics) 0, classical trajectories emerge from the lin- 
ear Schrodinger equation and the two types of quantum 
statistics have been linked to different classical behaviors: 
Poisson to regular motion, while Gaussian Ensembles to 
chaotic motion. It is natural to ask whether these find- 
ings for the linear Schrodinger equation can be general- 
ized to other types of wave equations . The Bogoliubov 
equation Q obtained from the linearization about the 
ground state of the Gross-Pitaveskii (G-P) equation has 
a purely real spectrum, and there is also a classical limit 
in the sense of geometric optics. It therefore makes sense 
and will be very interesting to explore the relationship 
between the Bogoliubov level statistics and regularity of 
the corresponding classical trajectories. 

There is, however, an important difference between the 
two types of equations; while the Schrodinger equation 
is Hermitian, the Bogoliubov equation is non-Hermitian 
and its statistics can not readily be predicted by standard 
random matrix theory. In fact, the Bogoliubov equation 
belongs to the category of symplectic problems, describ- 
ing linearized motion about stationary states in nonlin- 
ear classical Hamiltonian systems. This can be easily 
understood by noting that the G-P equation does have 
a classical Hamiltonian structure (of infinite dimensions, 
though) and that the Bogoliubov equation describes 
excitations about a stationary solution of the G-P equa- 



tion. The non-Hermiticity of the Bogoliubov equation 
makes it allowable to have complex eigenvalues in gen- 
eral, which signifies instability of the stationary solution. 
This will not happen about the ground state (lowest en- 
ergy state) which is always stable. Therefore, our inves- 
tigation of the Bogoliubov problem should shed light on 
the behaviors of motions around stable stationary states 
in extensive classical Hamiltonian systems. 

In this Letter, we investigate the level statistics of Bo- 
goliubov elementary excitation in separable circular as 
well as nonseparable stadium billiards. These are the ex- 
citations of a system of interacting particles in contrast 
with the modes of non-interacting particles described by 
the linear Schrodinger equation 8, 9]. The classical tra- 
jectories of Bogoliubov waves are found to be regular in 
circular billiards and chaotic in stadium billiards. By 
mapping the non-Hermitian Bogoliubov operator to a 
real symmetric matrix, we find the mean field interac- 
tions in the condensate do not change the level statistics 
of Bogoliubov excitations. This surprising result is tested 
numerically by using a phase shift method for calculat- 
ing the excitation energy. In the regime of strong inter- 
action and low excitation energy (phonon), we map the 
Bogoliubov equation to an equivalent Schrodinger equa- 
tion with Neumann boundary condition and show that 
the statistics of Bogoliubov levels are the same as that 
for the Schrodinger equation, although interactions in the 
condensate do change the average number of levels up to 
a certain energy. 

Gonsider condensed atoms confined in a quarter- 
stadium shaped trap of area A (with length of the top 
straight side L, radius of the semicircle R) and height d, 
where d << R so that lateral motion is negligible and 
the system is essentially two dimensional ^3 • With only 
a quarter of a stadium, one is restricted to a single sym- 
metry class of the full problem . The dynamics of the 
BEG are described by the G-P equation 

= V'V^ + 5^I^N, (1) 
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Figure 1: Plot of phase shift 5 versus k/^/g'. 

where g = 2\/2a/d is the scaled strength of nonlinear 
interaction, N is the number of atoms, a is the s-wave 
scattering length. The ground state of BEC can be writ- 
ten as "0 = ipo {f) exp (—i fit), where fi is the chemical 
potential and ipo (r) can be taken as real. The length 
and the energy are measured in units of 6 = and 
jm}?' respectively, so that the scaled area of billiards is 
A = 7r/4. The dynamics of the elementary excitations are 
obtained by linearizing G-P equation about the ground 
state and their energy spectrum is described by the 
time-independent Bogoliubov equation 



= E 



Hi H2 
H2 Hi 
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where is the Pauli matrix. Hi = +2gNil;Q — fj,, 

H2 — gNipQ, E is the Bogoliubov excitation energy, and 
(m, v) is the eigenwavefunction of linear operator C. The 
ground state wavefunction ipo (excitation {u,v)) has the 
normalization J^ipQdxdy = 1 {j ^{v? — v'^)dxdy = 1) and 
satisfies the boundary condition V'oaA — ^ {{u, v)^^ — 0). 

Classical trajectories or rays arise from the Bogoli- 
ubov equation in geometric optics approximation 
Assume a trial Bogoliubov wave of the form {u, v) = 
(a, (3) e^^ and consider a slowly- varying medium (Va, V/3 
small) and a slowly- varying velocity (V^S* small) approxi- 
mation. We obtain the Eikonal equation |VS'| = with 

p = \jiJL + {gNipof - 2gN^^. The classical tra- 

jectories of the Bogoliubov waves are still governed by 
the ray equation ^ (pw) = Vp, where w is the direc- 
tion of the trajectory and w is the arc-length coordinate 
along the trajectory. Interestingly, the classical trajecto- 
ries are simply straight lines for non-interacting as well 
as interacting uniform gases because p = \fEk, the lo- 
cal momentum in both cases. In the regime of strong 
interaction where the ground state of BEC is nearly uni- 
form, the classical trajectories of Bogoliubov waves are 
straight lines and undergo elastic specular reflection law 
at the boundary of the billiard. Therefore we predict that 
the Bogoliubov level statistics are still Poisson in circular 
billiards and Gaussian Orthogonal Ensembles (GOE) in 
stadium billiards through quantum classical correspon- 
dence. 

This prediction is supported by a general argument 



based on mapping the non-Hermitian Bogoliubov oper- 
ator £ to a real symmetric matrix. The linear operator 
£ can be written as £ ~ <^zQ, where Q is a real sym- 
metric positive definite matrix because the ground state 
of BEC is thermodynamical stable . The positive def- 
initeness of Q yields the decomposition Q = T'^T (T is 
a real matrix with nonzero eigenvalues) and the Bogoli- 
ubov equation reduces to 



E T 



(3) 



Therefore Bogoliubov excitation energy is the eigenvalue 
of a real symmetric matrix Tcr^T^ and should have GOE 
distribution in stadium billiards for arbitrary interaction 
strength 3]. 

In the following, we report numerical test of this pre- 
diction by developing a phase shift method to calcu- 
late the Bogoliubov excitation energy. Notice that the 
condensate density is nearly uniform in the inte- 
rior of billiards that yields the planewave forms of the 
Bogoliubov excitations. Similar to the scattering wave 
method in linear quantum mechanics [l^ . the nonuni- 
form condensate density close to the boundary and the 
hard walls can be taken as a pseudopotential and the 
scattering by this pseudopotential only induces a phase 
shift of the excited planewave in the interior of the bil- 
liard. The phase shift may be determined by solving 
the one-dimensional G-P equation with an infinite wall 
at a; < 0. Far from the wall, the Bogoliubov equation 
has both planewave and exponential solutions for a cer- 
tain excitation energy. We numerically integrate |l3l | the 
one-dimensional Bogoliubov equation with two different 
initial conditions to eliminate the exponential terms and 
extract the planewave solutions. The phase shift S is 
obtained by comparing the numerical solution with the 
expected sin(A:a; -I- S) dependence. The result is shown in 
Fig. 1 as a function of kj^/^ , where g' = gNip^, and (/?o 
is the condensate wavefunction far from the wall. We see 
that the phase shift approaches 7r/2 in the limit of low 
energy excitation and strong interactions (phonon) and 
asymptotically approaches zero in the regime of high ex- 
citation energy and weak interaction (free particles). As 
expected, the phase shift is always zero for g' — Q. 

With the phase shift method, we can calculate the Bo- 
goliubov excitation energy in one-dimensional billiards. 
The reflected planewaves from the boundary walls a; = 
and L can be written as (pQ = C sin (kx + 6k) and (f>L = 
Fsin(— fc(x — L) + 6k), respectively. The continuum of 
the wavefunction and its derivative in the interior of the 
billiards require (po/p'o = 4>l/4''l that yields the quanti- 
zation condition 



kL + 26k 



(4) 



where n is an integer. Eq. (4) determines wavevector 
k and Bogoliubov excitation energy E = J ^ (j^ + 2g') 
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Table I: Comparison of Bogoliubov excitation energies in one 
dimensional billiards using both phase shift (Ei) and matrix 
diagonalization methods (-E2). gN = 1000. 



El 





105.8 
212.2 
320.1 
430.4 
543.6 
660.4 
781.5 



2E-12 
102.6 
213.6 
322.0 
432.4 
545.6 
662.5 
783.6 



El 



En 



907.4 
1038.7 
1175.9 
1319.6 
1470.0 
1627.0 
1791.7 
1964.1 



909.7 
1041.1 
1178.3 
1321.9 
1472.3 
1629.7 
1794.5 
1967.0 



Si 



Eo 



2144.4 
2333.0 
2529.7 
2735.3 
2949.0 
3171.9 
3403.6 
3644.1 



2147.4 
2336.0 
2532.8 
2738.2 
2952.3 
3175.1 
3406.9 
3647.6 



El 



Eo 



3894.3 
4153.1 
4421.4 
4699.1 
4968.1 
5282.2 
5588.4 
5903.7 



3897.5 
4156.6 
4424.9 
4702.6 
4989.6 
5286.0 
5591.9 
5907.3 



in one-dimensional billiards. For the non-interacting case 
((7 = 0), the phase shift ^ = and Eq. (4) reduces to fc = 
n7r/i, the quantization condition for a single particle in 
one-dimensional billiards. As k approaches zero {E 0), 
the phase shift 5k 7r/2, therefore fc = (i? = 0) is a 
solution of Eq. (4) that corresponds to the ground state 
of BEC (the uniform density in the interior of the billiard 
indicates the 7r/2 phase shift). 

To check the validity of the phase shift method, we also 
calculate the Bogoliubov excitation energies using tradi- 
tional matrix diagonalization method in which the linear 
operator C is represented as a matrix and the diagonal- 
ization process gives the excitation energy E. The results 
are compared with those from the phase shift method in 
Table 1. We see that the phase shift method gives ac- 
curate results for the Bogoliubov excitation energy. The 
wavefunction of the first three excited states is shown in 
Fig. 2. Clearly, the excitation wavefunction is described 
by the planewave sin {kx -\- 6) in the interior of the bil- 
liards and drops to zero at the boundary. 

In the regime of weak interaction and high excitation 



as 




Figure 2: Bogoliubov excitation wavefunctions in one dimen- 
sional billiards for gN = 1000. Solid lines are from the matrix 
diagonalization method and dotted lines are from phase shift 
method smikx -I- S). (a), (b) and (c) represent the first three 
excitation wavefunctions, respectively. 



energy (A;^/2 >> 2(7'), the excitations behave like free 
particles where the spectrum is well understood. Di- 
rect calculations of Bogoliubov excitation energies in two- 
dimensional billiards for arbitrary interaction strength 
and excitation energy are difficult for both phase shift 
and matrix diagonalization methods. However, we are 
more interested in the regime of strong interaction and 
low excitation energy (A;^/2 << 251', phonon), where the 
effect of interactions is essential. In this regime, the phase 
shift b is approximately 7r/2 as seen from Fig. 1, which 
means that the first derivative of the planewave should 
be zero at the boundary, instead of the zero wavefunc- 
tion for the single particle case. The Bogoliubov wave- 
function far from the wall can be written as {uk,Vk) = 
(Uk, Vk) , where {Uu^Vk) = \{x + X-\x~ X^^) , X = 



1/4 



Wj2+2g 

with Neumann boundary condition 



satisfies the Schrodinger equation 



dn 



0, 



(5) 



71 is the normal direction of the billiard wall. Eq. (5) 
has an analytical solution (p = BJ2m (kr) cos {2m9) in 
quarter-circular billiards, where (kr) is the Bessel 
function, m is an integer, and B is a normalization 
constant. The boundary condition is satisfied by re- 
quiring Jjm (kR) — that determines wavevector k 
and Bogoliubov excitation energy. No analytical solu- 
tion is available for stadium billiards and we use the 
ansatz of the superposition of planewaves </)(a;, y) = 
X]j=i cos(kjxx) cos{kjyy) (planewave decomposition 
method , where kj^ — kcos{9j),kjy = fcsin(6'j), M 
is the number of planewaves, 9j — 2j7r/M, i.e. the direc- 
tion angles of the wave vectors are chosen equidistantly. 
The ansatz solves Eq. (5) inside the billiards and the 
boundary condition determines the possible wavevector 
k and Bogoliubov excitation energy. 

The average number of energy levels (N{E)) up to E 
should satisfy a Weyl-like type formula [l5j 



{N{E)) 



1 

in 



AE' + v^TEy + c 



(6) 
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Figure 3; Plots of the average number of energy levels (A'^ (E)) 
up to energy E in quarter-circular (a) and quarter-stadium 
(b) billiards, g' = gNipQ — 25000. Solid and dotted lines are 
from numerical results and Eq. (6), respectively. 
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Figure 4: The distributions of the spacing of neighbor- 
ing Bogoliubov energy levels ( k < 100 and the lower 
50 ones are omitted) in phonon regime, (a) Circular bil- 
liards. Dashed line: Poisson distribution Po{s) = exp(— s). 
(b) Stadium billiards. Dashed line: Brody distribution 
P (s) = (z^ -I- 1) a,^s'' exp (— a^s^^^) from the prediction of 

GOE, where a^, — [F (^j^pf )] ^ and F is Gamma function. 
The fitting Brody parameter f = 0.76 is smaller than ex- 
pected 0.953 due to the "bouncing ball states" "P] and finite 
levels effects [T^. 



where E 



fc2 




A and V are 



the area and the perimeter of the biUiard, and C is a 
constant related to the geometry and topology of the 
billiard boundary. Eq. (6) is only valid in the semi- 

that yields the Bo- 

In the limit 

1), we have 



classical limit E — 

E 
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goliubov energy ^ ^ , , -j^ 
of large interaction constant (Ejg' 

47r V .9' 2(g') 
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Cjwith the 



condition ^ < -p= < 

We have computed the Bogoliubov excitation ener- 
gies in circular billiards {R — 1,L — 0) using the an- 
alytical solution and in stadium billiards {R = L — 
ll + 4/7r)) using planewave decomposition method. 
In Fig. 3, we plot {N [E)) in both circular and stadium 
billiards. The agreement between the numerical results 
and the corresponding Weyl-like formula (Eq. (6)) is sat- 
isfactory. We unfold the spectrum formed by the energies 
Em i.e., we evaluate Eq. (6) for each En in order to ob- 
tain the new energies En — {N{En))- Note that the 
integer part of En is about n and, as a result, the cor- 
responding mean level spacing is characterized through 

(s) = {Eu+i - I {N) « 1. The resulting level 
spacing distributions P{s) are shown in Fig. 4. Clearly, 
the statistics of Bogoliubov excitation energy levels spac- 
ing are still Poisson in circular billiards and GOE in sta- 
dium billiards. 

Experimentally, the system can be realized by confin- 
ing BEC in two-dimensional optical billiards [Fg • Atoms 
could be trapped in a one dimensional optical lattice that 
is in the vertical direction to counteract gravity. Trans- 
verse motion could be confined to a closed pattern of light 
that is tuned to the blue side of atomic resonance, form- 



ing a repulsive barrier for the atoms. This pattern can be 
created by rapid scanning of a beam as was demonstrated 
for ultracold atoms. However, for the case of a BEC it 
would be better to created a static patterns, which can be 
accomplished using a liquidcrystal spatial light modula- 
tor [13 • The interaction strength may be adjusted by the 
confinement in the vertical direction, or by the number of 
atoms in each node of the standing wave. The Bogoliubov 
excitation energy may be measured using Raman tran- 
sition between two hyperfine ground states, denoted |1) 
and 1 2) respectively. The condensate would be formed in 
state Two co-propagating Raman beams would drive 
a transition to state |2) and the number of atoms in that 
state would be measured as a function of the frequency 
difference in the beams. The coupling efficiency to the 
excited states must depend upon the spatial profiles of 
the beams, requiring more detailed analysis (19l] . 
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